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Abstract 

In this and subsequent paper jl| we develop a recursive approach for calculating the short-time 
expansion of the propagator for a general quantum system in a time-dependent potential to orders 
that have not yet been accessible before. To this end the propagator is expressed in terms of a 
discretized effective potential, for which we derive and analytically solve a set of efficient recursion 
relations. Such a discretized effective potential can be used to substantially speed up numerical 
Monte Carlo simulations for path integrals, or to set up various analytic approximation techniques 
to study properties of quantum systems in time-dependent potentials. The analytically derived 
results are numerically verified by treating several simple models. 
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I. INTRODUCTION 



Studying quantum systems in time-dependent potentials represents a fundamental prob- 
lem which emerges in many areas of physics. Even if the Hamiltonian of the system itself 
does not explicitly depend on time, such situations naturally occur in the presence of time- 
dependent external fields. Another important example are fast-rotating systems, where the 
rotation introduces an explicit time dependence of the respective variables in the co-rotating 
frame. This includes modern experiments with rotating ultracold atoms and Bose-Einstein 



2|-l9|, studies of vortices 
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as well as non- stationary optical lattices 



12 



condensates 

Q. In the latter case counter-propagating laser beams are not perfectly modulated, so 
they produce an effectively moving optical lattice. An explicit treatment of time- dependent 
problems is also often required in nuclear and molecular physics 15|, and in studies of en- 
tanglement phenomena 16l4l8l|. Note that some systems can also be treated effectively with 
a time- dependent potential, which actually represents higher nonlinear terms in the wave 
function. A prominent example is the Gross-Pitaevskii equation, where the nonlinear term 
is treated within the split-step Crank- Nicolson method as a time-dependent potential which 
is updated after each time step |l9 |. 

A time- dependent formalism has been developed in the framework of various theoreti- 
cal approaches, based on the time-dependent Schrodinger equation. The well-known time- 
dependent semiclassical approximation is formulated by the Van Vleck propagator, and 
provides wide-range of important results, including generalizations to the time-dependent 
variational perturbation theory. Time-dependent variational principle, formulated by Dirac, 



represents basis for the time-dependent self-consistent fields method [20|. Another general 
method is time-dependent perturbation theory, which is also applied to derive a number of 
other approaches, including time- dependent scattering theory, linear response theory (Kubo 
formula) and Fermi's golden rule. Path integral formalism 2l|] also provides general frame- 
work to study dynamics of time-dependent quantum systems. In addition to these general 
methods, time- dependent variants of many highly specialized methods are also developed. 



As important exam 



Group (DMRG) 



es, let us mention time-dependent Density Matrix Renormalization 



231 ] ■ the Density 



Matrix Functional Theory (DMFT) [26 1. 



unctional Theory (DFT) 



24 



25| . and the Density 



In numerical approaches for time-dependent systems, partial differential equation repre- 
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sentation of the time-dependent Schrodinger equation can be solved using the discretization 
of the time domain through the finite difference method. In this case, the specific techniques 
include both explicit and implicit schemes, with the commonly used Crank-Nicolson scheme 

1 as the split- 



[27|. Other major numerical approaches include molecular dynamics, as we 
operator method, which is based on a space-grid calculation both without 



Monte Carlo 
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ail- 



28 



29| and with 



This method is mainly used in its second-order variant in the time of 



propagation for time-dependent potentials, where no change is required in comparison with 



the time- independent case. 



derived, including fourth 
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however, also higher order split-operator schemes have been 



35| and higher-order expansions 36l-l39[|. We also mention im- 
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42|, where discretization errors 



portant counterexamples, like world-line Monte Carlo 

i 

are completely removed [43!] , or open systems where the main contribution to the systematic 
error is due to retardation effects which cannot be dressed into the form of simple potentials 



44|. 
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In this and subsequent paper we extend the earlier established approach in Refs. 
48j for obtaining a high-order short-time expansion of transition amplitudes of time- 
independent potentials to the important time- dependent case. This development allows 
a high-precision calculation of transition amplitudes, which is necessary for extracting 
properties of various quantum systems, such as partition functions. Note that individual 
transition amplitudes can be accurately calculated using the lower-order effective actions at 
the expense of increasing the number of Monte Carlo time-steps, which would just increase 
the needed computation time. Although the presented approach can be used for improving 
the efficiency for calculating such transition amplitudes, it is mainly developed for applica- 
tions which require a large number of accurate transition amplitudes for further numerical 
calculations. Such situations occur, for instance, for determining partition functions 47. 14^ 



or for obtaining energ y sp ectra with the method of diagonalizing the space-discretized evolu- 



tion operator matrix 



50 



51| . In order to avoid an accumulation of numerical errors, so that 



such calculations can be performed with the required accuracy, the transition amplitudes 
have to be known more precisely. 

Our approach has already been successively used to study global and local properties 
of rotating Bose- Einstein condensates j52[, where the availability of accurately calculated 
thermal states has allowed to calculate precisely the condensation temperature, the ground- 
state occupancy, and time-of-flight absorption pictures even in the delicate regime of a critical 



and an overcritical rotation. If one needs a large numbers of accurate energy eigenvalues 
and eigenstates, as in this case, higher-order effective actions become useful. The precise 
knowledge of the short-time expansion of the propagator might also be used to improve 



28 



where 



numerical studies of time evolution, especially in systems exhibiting nonlinearities 
a time-dependent formalism is required. In addition to this, the presented approach can offer 
a significant improvement in studying dynamics of quantum systems within the Path Integral 



53| might be used to address the 



Monte Carlo calculations. The multilevel blocking method 
inherent dynamical sign problem, while the required number of Trotter time slices may be 
significantly decreased, due to high order of convergence of short-time propagators derived 
here. 



Thus, in spirit of Symanzik's improved action programme in quantum field theory [5j,|55|, 
we will introduce and analytically derive effective actions for time- dependent potentials 
which substantially speed up the numerical calculation of transition amplitudes and other 
quantities for such quantum systems. This approach does not suffer from technical problems 



related to backwards diffusion in time, as observed and discussed by Sheng 56(| and Suzuki 



571]. In the split-operator method one has to resort to multiproduct expansions 58| in order 
to resolve this problem. In the higher-order effective action approach for time- dependent 
potentials, which is presented here, the convergence of transition amplitudes is always guar- 



anteed in a natural way, as was shown conclusively for the time-independent case in Ref. 59 |. 
While this approach has a number of benefits (systematic higher-order expansion, simple 
implementation of symbolic calculation of effective actions using recursive relations, un- 
conditional convergence of transition amplitudes in the imaginary-time formalism, simple 
numerical implementation of the algorithm), it also has some disadvantages when compared 
to other short-time approaches. Most notably, higher-order effective action in this method 
contain higher spatial and time derivatives of the potential, and thus require the potential 
to be smooth function. Another possible disadvantage is increasing computational complex- 
ity of higher-order effective actions, which significantly grows with the order, and has to 
be appropriately chosen so as to minimize the computing time for the desired precision of 
numerical calculations. 

The outline of this paper is as follows. In Se c . ITIl w e briefly review the underlying 
path-integral formalism of quantum mechanics 
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6o|-62 



62| for systems in time-dependent 
potentials in order to fix our notation. Afterwards, Sec. [TTll presents a lowest-order path- 



integral calculation of the short-time transition amplitude for time-dependent potentials 
in the single-particle one-dimensional case. Based on these lowest-order results and using 
Schrodinger equations for transition amplitudes of systems in time-dependent potentials de- 
rived in Sec. IIV[ we develop in Sec. |V] the systematic efficient recursive approach for one 
dimensional quantum systems. The analytically derived results are verified numerically with 
the help of several illustrative one-dimensional models in Sec. IVII 

The approach developed here to obtain transition amplitudes for time-dependent poten- 
tials is generalized in the second part of this paper \l\ to quantum systems with many degrees 
of freedom. In the next paper we also show how the developed imaginary-time formalism 
is transformed into a real-time one, and demonstrate its applicability by treating several 
models. 

II. PATH-INTEGRAL FORMALISM FOR SYSTEMS WITH TIME-DEPENDENT 
POTENTIALS 

We will consider a non-relativistic quantum mult i- component system in d spatial dimen- 
sions with a Hamilton operator which consists of the usual kinetic term and a time-dependent 
potential: 



Here P stands for the number of particles, the P x d dimensional vectors q and p describe 
positions and momenta of all particles, and the parenthetic subscript (i) denotes the cor- 
responding quantity for particle i. Although in this paper we only consider single-particle 
systems, here we introduce many-body notation and use it in introductory parts of Sections 
mil and HV ] as a preparation for the second part of this paper where extensions of the 
presented approach to many-body systems and to the real-time formalism are developed. 

Note that the potential V(q, t) of the system is allowed to depend on the positions 
of all particles, and, therefore, contains implicitly all types of interactions. In practical 
applications the potential V usually contains, apart from the external potential, also two- and 
three-body interactions, but further many-body interactions can, in principle, be included 
as well. 




p 



2 



(2.1) 



i=l 
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The central object for studying the dynamics of such a quantum system within the path- 
integral formulation is the transition amplitude 

A(a, ta] b, h) = (b, h\U{ta ^ tb) |a, 4) . (2.2) 

Here the vectors a and b describe the positions of all particles at the initial and final time 
ta and tb, \si,ta) and \h,tb) denote the corresponding Hilbert-space states of the system, and 



U{ta^U)=fexp}^-^j^ dtH{p,ci,t)j (2.3) 
represents the evolution operator of the system describing its propagation from ta to t^. Here 



we assumed ta < tf, and introduced the standard time-ordering operator 
. . . , f d(t)d(t'), if t > t', 

T{0{t)0{t')} = <^ ^ ^ ^ ^' . (2.4) 

I d{t')d{t), otherwise. 

The starting point in setting up the path-integral formalism is the completeness relation 



A(a, ta] b, Q = J dqi-- J dq^N-i ^(a, ta] qi, ti) v4(qi, ti; q2, ^2) ■ ■ ■ ^(qAf-i, tN-i] b, tb) , 

(2.5) 

where e = (tb — ta)/N denotes the time-slice width, tn = ta + ne are discrete time steps, 
and the P x d dimensional vectors qi, . . . , qAr_i describe positions of all particles at a given 
discrete time step which is specified by the non-parenthetic index. To leading order in e, 
the short-time transition amplitude reads 

A{qn,tn]qn+i,tn+i) ^ cxp | ^g^^^ (q„, q„+i, t„+i) | , (2.6) 

where the naive discretized action S^^^ is usually expressed as 



(2.7) 



Here we introduced the discretized velocity Sn = qn+i — qn, and rescaled the coordinates 
so that the mass of all particles is equal to unity. The potential V is evaluated at the 
mid-point coordinate x„ = (q„ + q„_|_i)/2 and at the mid-point time r„ = (t„ + tn+i)/2. 
Eq. (12. 7p is correct to order e ~ and, therefore, after substitution to Eq. (12.61) . leads 
to errors of the order 0(1/A^^) for discretized short-time transition amplitudes. Note that 



the normalization factor ~ in Eq. f l2.6p does not affect the A^-scahng of errors, 

since short-time amphtudes will be inserted into the completeness relation (12.51) . where 
this normalization factor will be added to the normalization of the long-time transition 
amplitude. However, due to the fact that Eq. (12. 5p contains a product of short-time 
transition amplitudes, the deviation of the obtained discrete transition amplitude from the 
corresponding continuum result will be of the order N ■ 0{1/N'^) = 0{1/N). For this reason 
the naive discretized action is designated by S^^\ 

In the limit — i- oo we recover the continuous transition amplitude, which leads to the 
formal coordinate-space path-integral expression 

/.q(tb)=b 
Jq,{ta)=a 

where the integration is defined over all possible trajectories q(t) through the discretization 
process described above. In this equation, the action S for a given trajectory q(t) is defined 
as usual: 



A(a, ta, b, h) = I Vci{t) exp { ^S[q\ } , (2.8) 



1 



S[q]= rft |-q2(t)-r(q(t),t)| . (2.9) 
The outlined derivation represents the basis for the path-integral formulation of quantum 



mechanics 
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60|-|62| as well as for the numerical calculation of path integrals. The described 



discretization procedure is most straightforwardly numerically implemented by the Path- 



Integral Monte Carlo approach 



63| 



Note that the 0{1/N) convergence can be also achieved with other choices of space-time 
points at which we evaluate the potential in Eq. (12. 7p . For example, the left or the right 
prescription x„ = q„, r„ = t„ or x„ = q„+i, r„+i = tn+i is often used for the spatio-temporal 
argument of the potential. These different choices do neither affect the numerical calculation 
nor the analytical derivation, and different prescriptions can even be translated into each 
other. However, it turns out that the mid-point prescription always yields the simplest 
analytic results, and, therefore, we will use it throughout the present paper. 

In the following we will switch to the imaginary-time formalism, which is widely used 



in numerical simulations [63|], since it mitigates problems which are related to the oscil- 
latory nature of the integrand in the real-time approach. To obtain the real-time from 
the imaginary-time amplitudes, one would have to perform an inverse Wick-rotation, i.e. a 
suitable analytic continuation of the numerical results. This might be difficult due to the 



inherent instability of this procedure with respect to statistical noise, which is always present 
in numerical calculations. However, using the imaginary time is justified by the fact that 
all current applications of this approach are related to quantum statistical physics, which 
is naturally set up in imaginary time, with the inverse temperature /3 = l/Zc^T playing the 
role of the imaginary time. Also, energy spectra and energy eigenfunctions can be efficiently 
calculated in this formalism 
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51| . Note that the derived analytic expressions for higher- 



order propagators can be formally transformed from the imaginary- to the real-time axis. 
In the next paper we will demonstrate by treating several simple examples that such 
analytic expressions can be successfully used for calculating the real-time evolution in the 
context of the space-discretized approach jsol. Isi]. 

After Wick rotation to the imaginary time, the transition amplitude is expressed as 

A(a,4;b,4) = / I)q(t) e-TT^^^l'^l , 

J(l{ta)=a 

where the action is replaced by its imaginary-time counterpart, the Euclidean action 



(2.10) 



ta 



dt \]^^\t) + V{ci{t),t) 



(2.11) 



which represents the energy of the system. To simplify the notation, we will drop the 
subscript E from now on. If we consider in Eq. (12.101) only diagonal amplitudes, i.e. a = b, 
and integrate over a, we obtain the path-integral expression for the partition function 

-hp 



Z{/3) = Tr |Texp 
by setting ta = and tb - 



1 

h(3. 



dtHit) 



(2.12) 



III. PATH-INTEGRAL CALCULATION OF THE PROPAGATOR 

In this paper we develop a method for calculating short-time transition amplitudes for 
time-dependent potentials to high orders in the propagation time. To this end we follow the 
approach of Ref. 48| and generalize the level p = 1 discretized transition amplitude from 



Eq. (12. 6p in such a way that the exact transition amplitude is written as 

^(^'^-b,^^)= (^;^^"'*^"'"^^- (3-1) 



Here we use the convention h = 1, the ideal discretized effective action reads 59|, |64| 



S*{^,d-e,r) = ^ + eW{^,d-e,T), (3.2) 

and W represents the ideal effective potential, which ensures the exactness of the above 
expression. The latter depends not only on the coordinate mid-point x = (a + b)/2, the 
discretized velocity (5 = b — a, and the time interval e = — as already introduced in the 
case for time-independent potentials, but also on the time mid-point r = {ta + tb)/2, due to 
the explicit time dependence of the potential. 

We will analytically derive a systematic short-time expansion of the effective potential W, 
which provides an improved convergence for numerically calculating transition amplitudes 
and partition functions as well as other properties of quantum systems with time-dependent 
potentials. The expansion of the effective potential W to order e^~^ yields the effective 
action correct to order e^, with errors proportional to e^^^. Due to the normalization factor 
in the expression (13. ip . the total e- convergence of the amplitude is given by 

Ap(a, ta, b, tfe) = /l(a, 4; b, t^) + 0(^^+1-^/2) ^3 3) 

This e-scaling of errors is valid if we are interested in calculating short-time transition 
amplitudes, which is the main objective of this paper. However, if we use such short- 
time transition amplitudes to calculate long-time transition amplitudes through the time- 
discretization procedure (12. 5p . due to successive integrals the normalization factors will be 
again added, and we will get total errors of the order ■ 0{1/N'p^^) = 0{1/N^). Although 
the corresponding effective actions and resulting discretized short-time transition amplitudes 
are designated by the index p according to their A^-scaling behavior, we stress that the scaling 
with respect to the short propagation time is always given by Eq. (13.31) . 

Before we embark upon developing a systematic recursive approach for analytically cal- 
culating higher-order effective actions, we first have to study the general structure of the 
effective potential, which turns out to be more complex than in the case of time-independent 
quantum systems. In order to do so, we calculate the short-time expansion of the effective 



potential by using an ab initio approach introduced in Ref. [47j. To simplify the calculation, 
we will restrict ourselves in the present section to the single-particle one-dimensional case. 
Based on these results we will develop in Sec. |V] the systematic recursive approach for such 
simple quantum systems, and then extend and generahze this procedure to systems with 
many degrees of freedom. 



Following Ref. j47|, we start with changing the variables via q{t) = ^{t) +y{t), where ^(t) 
is some chosen reference trajectory satisfying the same boundary conditions as the path q{t), 
i.e. ^(ta) = O', ^{tb) = b. This implies that the new variable y{t) vanishes at the boundaries, 
i.e. y{ta) = y{tb) = 0. We also introduce a new time variable s by t = r + s, in which the 
boundaries are defined by Sa = — £^/2, = s/2. In the new variables, the kinetic energy 
functional has the form 

where a dot represents a derivative with respect to the time s. With this the transition 
amplitude reads 

A(a,t„;6,t,)=e-^K1W / e-/i/2'^4l^^(^H^efeW.«)} , (3.5) 

Jy{-e/2)=0 

where the quantity is defined as 

s) = Vi^s) + y{s), T + s)- Vi^s), t) - y{s)i{s) , (3.6) 

and the Euclidean action 5'[.^](r) for the reference trajectory C,{s) is defined by an expression 
similar to Eq. (12. lip , but with the time argument of the potential being fixed now by the 
mid-point r: 

Smr) = f'\s ^-i\s) + V{as\r)^ ■ (3-7) 
Thus, the transition amplitude reads 

\/2'Ke \ / 
where the path-integral expectation value is defined with respect to the free-particle action: 

ry'-^/^)=^ re/2 ,12,, 

{•) = V2^ Vy{s) • e-^-^/^"'-^"^ '^'K (3.9) 

A{-e/2)=0 

The transition amplitude ( 13. 8 p can then be obtained through a standard calculation of the 
free-particle expectation value by using the Taylor expansion 

re/2 



1 



e/2 1-6/2 



+ -/ ds I ds'{U^{y{s),s)U^{y{s'),s')) + ... .{^.IQ) 



e/2 J -e/2 
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As we see from Eq. ( 13.60 . the quantity has the simplest form if we choose the reference 
trajectory ^ in such a way that its second derivative with respect to the time s vanishes. 
The natural choice is thus the linear trajectory 



^{s) = x + s6/e, 

centered around the mid-point x = {a + b)/2. In order to calculate the expectation values 
in Eq. fl3.10p . we further have to expand 

UM^), s) = Vi^s) + y{s), T + s)- Vi^s), t) (3.11) 

around this reference trajectory according to 

oo 1 (■ "I 

UMs),s)= ZJZr, V^-\as),r)y\s)s-, (3.12) 

Tl. Til. 

n.m — \.) 
n-|-m>0 

where, in order to simplify the notation, (m) denotes the order of the partial derivative 
with respect to the time s, and (n) denotes the order of the partial derivative with re- 
spect to the spatial coordinate. For a free-particle theory, expectation values of the type 
(y"^(si) y"'^{s2) . . . ) can be easily calculated using the standard generating functional ap- 
proach. The expectation value {y{s)) vanishes due to the symmetry, while the correlator 
A(s, s') = {y{s) ?/(s')) is given by the expression 

A(., .') = (I - ^) (I + s') + (. o .') , (3.13) 



and higher expectation values can be found in Ref. [4j 

In order to obtain an expansion of the transition amplitude in the propagation time e, one 
has to consider Eq. f l3.10p and to take into account the powers of e of all terms to identify 
the relevant terms at the desired level p. To illustrate this, let us look at the term linear in 

r^\s{UM^),s))= f'\s -^y^-\as\r){y%s))s-, (3.14) 

J-e/2 J-e/2 ^' ^' 

n + m > 

where the integration over s yields terms proportional to e™'^^ times the contribution of the 
expectation value term. For example, from Eq. fl3.13p we see that the n = 2 expectation 
value would yield an additional e power, amounting to a total dependence of the term 

11 



corresponding to a given m and n = 2. Since we are interested to calculate all terms to a 
given order in e, we have to carefully inspect all possible terms and to select the appropriate 
ones at a given level p. As an example we write down all terms which contribute to the 
expectation value f l3.10p to order 0{e'^): 



^^-rS^,.sU,iyis),s)^ = 1 - r ds {lv"{as),r){y\s)) + ^^V^'\as),r){y\s)) 



e/2 

-e/2 

+ Vi^s), t)s + ^V^'Kas), r){y'{s))s + ^Vi^s), r)s' 
+ \ r r ^'(^(^)' ^)^'(^(^'), r) {y{s)y{s')) + 0{e') . (3.15) 

^ J -e/2 J-s/2 

In the above expression, the correlators = A(s, s) and {y{s)y{s')) = A(s, s') are given 

by Eq. f lS.lSp . and the expectation value can be directly determined using either the 

generating functional method or the Wick rule, yielding 3A^(s, s). In order to be able to 
calculate the remaining integrals over s and s', we need to expand also the potential V and 
its derivatives with respect to the first argument ^(s) = x + s6/e around the mid-point 
X. The required number of terms in this expansion, contributing to the desired order of 
e, is obtained by taking into account the diffusion relation 6"^ ~ e, which is valid for small 
propagation times e and has been demonstrated to yield a consistent expansion of short- 



time transition amplitudes 48|. The expansion of the potential into power series in s6/e 
gives an additional polynomial s-dependence, which finally allows an analytic calculation of 
all integrals in Eq. f l3.15p . When this is done, we obtain the following expression for the 
expectation value flS.lOp 

\ / 12 12 240 24 24 480 

_ v(3)M_y(3)^ + o(e^), (3.16) 
240 480 ^ ^' ^ ^ 

where the potential V as well as its spatial and temporal derivatives are evaluated at the 
mid-point x, r, i.e. V = V{x, r). Note that we have retained in Eq. f l3.16p only those terms 
whose order is less than e^. Taking into account the diffusion relation 5^ ~ e, the terms 
proportional to and are considered to be of the same order. 

Combining Eqs. (13. ip and (13. 8p . we see that the ideal effective action can be expressed 
by its short-time expression (13. 7p and the expectation value (I3.10p according to 

S*ix,5;e,T) = S[e](r) - log (e"^-/^'^^^^'^^^)'^)^ . (3.17) 
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After expanding the potential V around the mid-point x in the action 5'[^](r) in Eq. (13.71) . 
we obtain its short-time expansion: 

Sim = Y,+y^ + + V^"^^ + 0{e') . (3.18) 

This allows us to calculate the short-time expansion of the effective potential W from 
Eqs. (13. 2 p and (I3.16p - (l3.18p . For example, up to level p = 3 we get 



5' 



12 24 12 240 24 24 480 1920 

(3.19) 

Numerically, such a result allows to a speed up the calculation of transition amplitudes, since 
the errors can be substantially reduced using analytic expressions for higher level effective 
actions. 

Compared to our previous results for effective actions of time-independent potentials V 



in Ref. 481], we see that new terms appear which contain time derivatives of the potential. In 
particular, we observe the emergence of terms with odd powers of the discretized velocity 5, 
which was previously not the case. In fact, for time-independent potentials we have shown 
in Ref. 48|| that the effective potential is symmetric in 5, which leaves only even powers of 5 
in its short-time expansion. Here, however, also odd powers of 5 survive due to the explicit 
time dependence of the potential. We also recognize that all new terms are proportional to 
time derivatives of the potential and vanish in the time-independent case, thus reducing the 



48| 



effective action to the previous expressions in Ref. 

Therefore, the correct systematic short-time expansion of the effective potential turns out 
to have the form 



cxj lu I / r\2k / J\ \ 2fc+l 



m=0 A;=0 



(3.20) 



Here 5/2 is used as the expansion parameter in order to have expansion coefficients c which 
are defined consistently with Ref. |48]. Such an expansion allows that the level p effective 
action is written as the sum of terms corresponding to m = 0, 1, . . . ,p — 1. Note that for 
m = p — 1 we need to take into account only the even-power terms Cp_i^fc(x, r) e'^^^^'^ [5 / . 
The odd-power terms Cp_i/2,fc(a^5 t) ^''^^^'^(5/2)^^"*"^ are proportional to e:^^^/^, i.e. they are 

13 



of higher order than the required ^ for level p effective action. For this reason, the correct 
expansion of the effective potential at level p is given by 

iy,(a;,5;5,r) = J] J]c^,.(a;,r)5--M- + J] c„+v2,fc(x, r) e^'M - , (3.21) 

m=0 fc=0 ^ m=0 fc=0 ^ ^ 

and it provides the convergence of transition amplitudes according to Eq. (13. 3p . 



IV. FORWARD AND BACKWARD SCHRODINGER EQUATION 

In this and in the next section we will use the latter result (I3.20p to develop a systematic 
recursive approach for analytically calculating effective actions to high levels p for time- 
dependent potentials. To this end, in this section we rederive both the forward and the 
backward Schrodinger equation for the transition amplitude and then use them to derive 
corresponding differential equations for the effective potential W . Afterwards, in Sec. IVl we 
use Eq. (I3.20p to solve the derived equations for W and to obtain recursion relations for 
single-particle one- dimensional systems. 

The evolution operator for quantum system in a time-dependent potential is given by 
Eq. (12. 3p . or, in imaginary time. 



b 



U{ta^U)=Texp^- J dtH{p,q,t)j . (4.1) 
Thus, the evolution operator obeys the differential equation 

dt, U{ta ^ tf,) = -H{p, q, U) U{ta ^ h) , (4.2) 
and, similarly, 

dt, Uita ^ h) = Uita ^ h) H{p, q, ta) . (4.3) 

If we determine from Eq. (14. 2 p the matrix elements which correspond to the transition 
amplitude Eq. (12. 2p . we obtain the forward Schrodinger equation for the transition amplitude 

dt, A{a, ta, b, U) = -Hb A(a, 4; b, U) , (4.4) 

where Hb stands for the coordinate-space Hamilton operator Hb = H{—idh,h,tb), in which 
momentum and position operators are replaced by their coordinate-space representations at 
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b. To obtain the analogous equation for the derivative with respect to the initial time ta, 
we have to take into account that the imaginary-time transition amplitudes as well as its 
time derivative are real. With this Eq. fl4.3p yields at first 



Since we have in addition 



we finally get the backward Schrodinger equation for the transition amplitude 
dt, A(a, ta, b, tb) = Ha A{a, ta] b, U) , 



(4.5) 



(4.6) 



(4.7) 



where Ha = H{—ida, 8i,ta) is defined analogously as Hb- 

In the next step we change the original time variables ta and tb to the mid-point r and 
the propagation time e, which converts (14.41) and (14.71) to 



d, + -{Ha + Hb 



A{a,ta;h,tb) = 
dr + {Hb-Ha) A(a,t„;b,t6) = 



(4.8) 
(4.9) 



Subsequently, we also change the spatial variables a and b to the mid-point x and the 
discretized velocity x = 5/2, thus Eqs. (14. 8 p and (14. 9 p read then 



dr-\dd + v^ 



v4(x,x;£:,r) 
v4(x,x;£,r) 








(4.10) 
(4.11) 



where we have introduced V± = (x ± x, r ± |) as abbreviations, and the P x ci-dimensional 
Laplacians ■ = d"^ over coordinates x and d^ - = W' over coordinates x, as well 
as the mixed Laplacian ^x ■ = dd. Note that we will not use Eq. (14. lip in our further 
calculation of the short-time transition amplitude. This equation describes the dynamics 
caused only by the presence of the explicit time-dependence of the potential, as is indicated 
by the derivative with respect to the time mid-point r. 

If we now express the transition amplitude (13. ip by using the effective potential in (13. 2p 
and the new spatio-temporal variables, we get 

1 



A(x,x;e,r) 



(27r£)™/2 



g-fx2-eiy{x,x;e,T) 



(4.12) 
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Substituting this expression into the differential equation for the transition amphtude fl4.10p . 
we obtain a corresponding differential equation for the effective potential: 



W+^-dW+edeW-led^W-led'^W+le\dWf+le\dWf = {V++V-) . (4.13) 
8 8 8 8 2 



This equation is formally identical to the corresponding equation (29) from Ref. |48| . 
and in the limit of the time-independent potential we recover the previously derived result. 
Eq. 04.131) can now be used to develop a recursive approach for calculating the effective 
potential W by using the double power series fl3.20p in both e and x. 



V. RECURSIVE CALCULATION OF THE PROPAGATOR FOR ONE-DIMEN- 
SIONAL SYSTEMS 

In this section we consider a single-particle one-dimensional system in a given external 
time-dependent potential V, and fully develop recursive approach for calculation of the short- 
time expansion of the effective potential W. In this case, the right-hand side of Eq. ( I4.13p 
can be expanded in a double power series of the form 

k,m>0 

We see that m and k must be of the same parity in order to yield a non-zero contribution. 
Therefore, we introduce the quantity n(m, k), which is equal to one if m — k is even, and 
vanishes otherwise. This allows us to reorganize terms of the double power series (15. ip in a 
form which corresponds to the expansion of the effective potential in Eq. (I3.20p : 

^ ^ r U{m,k)e^-'x"' [l-n(m,fc)]5— %%k+i)\ .^ 

^^f^^\{2k)\{m-k)\2^-'' (2A; + l)!(m-A;)!2— ^ J • l^- J 

In order to solve Eq. f l4.13p with (15. 2p . we now substitute the effective potential W with its 
double expansion (I3.20p . By comparing terms with even and odd powers of x, we obtain 
algebraic equations which determine the coefficients Cm,k and Cm+i/2,fc- After a straight- 
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forward calculation we get 

{ra—k) 

Yiim k) V 

8(m + A: + 1) c^,, = § (3^), + (2A: + 2){2k + 1) c^,.+i + C,,, 

~ '^m-'-2,A:-r + '^i+l/2,r '^m-i-5/2,A:-r-l + (^^ — 2r + 2) Q^^ Cm-l-l,k-r+l 

l,r 

+ (2r + 1) {2k - 2r + 1) C/+i/2,r c„_,_3/2,fc-r} , (5.3) 

\l - Yiim k)] V 

8(m + + 2) c^+i/2,fe = 8 (2A: + l)!'(m -A:)!2"--fe + + ^^^^^ + c™+i/2,fc+i + 

- XI {^^A' C™_i_3/2,fc~r + Cj+l/2,r C^-i-2,fc-r + (2^ + 1) (2A; - 2r + 2) Q+i/2,r- Crn-i-l,fc-r+l 



+2r (2A; - 2r + 3) q,^ Cm-/-i/2,fc-r+i } • (5.4) 

In the above relations, sums over / and r are restricted by the non-negativity of in- 
dices for both even and odd-power coefficients c. It is immediately seen in the case of a 
time-independent potential that the recursion for even-power coefficients Cm,k in Eq. (15. 3p 

48|, while the odd-power recursion (15. 4p 



reduces to the previously derived recursion of Ref. 
automatically renders all coefficients Cm+i/2,k to be zero. 

In order to efficiently solve these algebraic equations at a given level p, we recall the 
corresponding hierarchical level p expansion (I3.2ip . In such an expansion, we look only for 
even-power coefficients Cm,k with < m < p — 1, and odd-power coefficients Cm+i/2,fc with 
< m < p — 2, while the index k is always restricted by < A; < m. To simplify the 
calculation, we formally set all the coefficients with higher values of indices to zero. As in 
case for a time-independent potential in Ref. [48|, we are then able to calculate explicitly 
diagonal even-order coefficient Cm,m for any m: 

y{2m) 



(2m + 1)! 

Correspondingly, we obtain for the diagonal odd-power coefficients the general result 



(5.5) 



Cm+l/2,m — . (5.6) 

Thus, the recursion relations (15. 3p and (15. 4p represent together with (15. 5p and (15. 6p a closed 
set of algebraic equations, which completely determines the coefficients c for the effective 
potential at a given level p. 
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To illustrate how the recursive procedure works in detail, we solve the single-particle 
recursive relations up to order p = 3 and demonstrate that we obtain the same result 
(I3.19P as in the previous section, where we used the path-integral approach to calculate the 
effective potential. For the level p = 3 we need coefficients c with the first index m in the 
range 0, . . . , p — 1, i.e. m = 0, 1, 2. For a given m < p— 1 we start with calculating even-power 
coefficients Cm,k by using the recursive relation (15.31) in the order k = m,m — 1, . . . ,0, where 
we take into account (15.51) . Then we proceed with the calculation of odd-power coefficients 
Cm+i/2,k by using Eq. (15. 4p for the order k = m — 1, . . . ,0 and by taking into account (15.61) . 
We also recall that odd-power coefficients need not be calculated for m = p — 1, i.e. in the 
last step of the recursive procedure for m = p — 1 we only calculate even-power coefficients. 
Following this procedure, we immediately from get (15. 5p for m = 

cofl = V. (5.7) 

The coefficient Ci/2,0 = follows trivial from (15. 6p . For m = 1 we get the following even- 
power coefficients, corresponding to the values k = 1 and k = 0: 

Ci,i = V, (5.8) 


Ci,o = ^<o + ^ci,i = ^\/". (5.9) 
For the only non-trivial m = 1 odd-power coefficient C3/2.0 we get 

C3/2,0 = \c,/2,l + ^c;'/2,0 + = ^V' . (5.10) 

Similarly, we find the m = 2 even-power coefficients by using Eqs. (15.30 and (15. 5p : 



^^. = ^y'". (5.11) 

C2 = — - —c'\ + —c'' +—C2l = —V+- — . (5.13) 

2,0 24 24 24 12 ' 24 240 24 ^ ' 

This result is already sufficient to construct the effective potential Wp=3. If we insert the 

calculated coefficients in the expansion (13.210 at the given level p = 3, we reobtain the same 

expression as in Eq. (I3.19p . Comparing the calculated even-power coefficients Cm,k with the 

previously calculated ones for the time-independent potential in Ref. [48[, we see that they 

coincide if we set all time derivatives of the potential to zero, as expected, and that all 

odd-power coefficients vanish. 
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We stress that the recursive approach is far more efficient than the path-integral cal- 
culation presented in the previous section. For example, if one wants to extend a level p 
calculation to a higher order p' = p + 1, this requires in the path- integral approach not 
only to take into account the next term in the expansion fl3.15p . but also each previously 
calculated expectation value term has to be redone to one order in e higher. The complexity 
of this algorithm prevents its efficient implementation. However, in the present recursive 
approach, all we need to do is to calculate one additional order of odd-power coefficients 
Cp_i/2,fc, which corresponds to m = p' — 2 = p — 1, and even-power coefficients Cp^k, which 
corresponds to m = p' — 1 = p. To do this, we just apply the recursive relations and use 
previously calculated lower-order coefficients c. For instance, in order to obtain level p' = 4, 
we proceed ffist with calculating m = p' — 2 = 2 odd-power coefficients C5/2,fe- The high- 
est coefficient 05/2,2 = is automatically equal to zero due to (15.61) . while for lower- level 
coefficients we get 

C5/2,l = ^C5/2,2 + ^4'/2,l + l/^'' = ^^V^'^ , (5.14) 

3 11 1 1 • 

C5/2,0 = Y^C5/2,1 + ^4/2,0 - Y^C;/2,oCo,0 " gCl/2,0Cl,l = ^^^^^^ ■ (5-15) 

To completely determine the p' = 4 effective action, we then calculate even-power m = 
p' — 1 = 3 coefficients cs^k, and obtain from Eqs. (15. 3 p and (15. 5 p 

C3,3 



5040 



y^*^) , (5.16) 

C3,2^^V(^ (5.17) 

C3 1 = ^^(6) + —V" - —V"^ - — rV(3), (5.18) 

' 3360 80 360 120 ^ ' 

6720 480 360 120 ^ ' 

The outlined procedure continues in the same way for higher levels p. We have automa- 



tized this procedure and implemented it in our SPEEDUP code j65|] using the Mathematica 
7.0 package 661] for symbolic calculus. Using this we have determined the effective action 
for a one-dimensional particle in a general time-dependent potential up to the level p = 20. 
Such calculation requires around 2 GB of memory and approximately 1.5 hours of CPU time 
on 2.33 GHz Intel Xeon E5345 processor. Although the effective actions grow in complex- 
ity with level p, the Schrodinger equation method for calculating the discrete-time effective 
actions turns out to be extremely efficient. The value of p = 20 far surpasses the previously 
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obtained best result known in literature of p = 6 [38[ , and is limited practically only by the 
sheer size of the expression for the effective action of a general theory at such a high level. 
The whole technique can be pushed even much further when working on specific potential 
classes as, for instance, polynomial potentials, where higher spatial and temporal derivatives 
have a simple form. However, if this is not the case, the obtained general expressions for 
the effective potential can be used for any given potential. 

In practical applications, the increasing complexity of expressions for higher-level effective 
actions leads to a corresponding increase in the necessary computation time. For levels 



P<10 



45 



461], the increase in the computation time is minimal compared to the obtained 



decrease in the error of numerically calculated transition amplitudes, thus it can be used to 
speed up the calculation at a given level of accuracy. However, for higher values of p the error 
may be more efficiently reduced by decreasing the time step, i.e. by using a larger number of 
Monte Carlo discretization steps. Since the complexity of effective actions strongly depends 
on the concrete form of the given potential, the optimal values of both the level p and the 
size of the time step have to be estimated from a series of scaled-down numerical simulations. 
For example, level p = 21 and a certain time step have turned out to be optimal for studying 
fast-rotating Bose-Einstein condensates in an anharmonic trapping potential 52 1. 

We will generalize the presented approach to many-body systems in the next paper 
which will also include generalization of the imaginary-time to the real-time formalism, 
suitable for study of dynamics of quantum systems. 



VI. NUMERICAL VERIFICATION 

In order to numerically verify the derived analytical results, we have studied several 
simple models with time-dependent potentials. We have used the modified version of the 



Monte Carlo SPEEDUP 



65l | code, which implements higher order effective actions in the C 



programming language. In order to be able to resolve decreasingly small errors associated 
with higher levels p of the effective potential, we decided to consider at first some exactly 
solvable potentials. 

In Ref. 



67( 1 it was shown that a generalized Duru-Kleinert transformation 
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69| 



allows to map the transition amplitude for a time-independent potential V{x) to the corre- 
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FIG. 1. Deviations of amplitudes as functions of propagation time e = tf, — ta for the harmonic 



oscillator (j6.2p with a; = 1, rescaled with the Grosche factor ({t) = \/t'^ + 1. (left) Deviations 
of amplitudes |^p(2.5, 0; 2.51, e) — j4(2.5, 0; 2.51, e)| as functions of e, calculated analytically for 
p = 2, 4, 6, 8, 10, 12, 14, 16, 18, 20 from top to bottom. The dashed lines are proportional to e^'^ 
and e'^^'^ and demonstrate the perfect scaling of the corresponding level p = 2 and p = 20 results, 
(right) Comparison of deviations |^p(l,0; 1.001, e) — vl(l,0; 1.001, e)| calculated using the correct 
level p = 1, 2, 3, 4 effective potentials (full lines, top to bottom) with the deviations obtained for the 



481 ] for the case of time-independent potential. 



same levels p of previously derived effective actions 
Deviations for different levels p of time-independent effective actions correspond to different point 
types. 



sponding transition amplitude for the time-dependent potential 
1 



VG{x,t) 



-V 



X 



(6.1) 



at) \m, 

where C^(^) = Co + Ci^ + C2^^ is a general quadratic polynomial in time. Fig. [T] presents 
the numerical results for the time- dependent harmonic oscillator potential, where the time 
dependence is introduced by using the Grosche rescaling factor ({t) = y/l + 1"^: 



Vc 



G,HO 



>,t) 



(6.2) 



2(l+t2)2 • 

In the left plot of Fig. [T] we see that the obtained discretized amplitudes converge to the 
continuum limit systematically faster and faster when we use higher level effective actions. 
This log-log plot demonstrates that the analytically derived law e^+^Z^ for deviations of 
single-particle discretized transition amplitudes in d = 1 from the continuum transition 
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amplitudes holds perfectly, which verifies numerically our analytical results. The deviations 
are calculated using the analytically known continuous transition amplitude 67|]. The graph 
on the right of Fig. [T] illustrates the importance of terms with time derivatives of the potential 
in higher-order effective actions. In this plot we show deviations of discretized transition 
amplitudes calculated using the effective actions for time-independent potentials, derived 
in Ref. ks 



and compare them with the deviations of discretized transition amplitudes 
calculated with the correct effective actions, derived here for the case of time-dependent 
potentials. As can be seen from the graph, time-independent effective actions do not improve 
results after level p = 2, due to the fact that terms, containing time derivatives of the 
potential, are not systematically eliminated from deviations when such effective actions are 
used. If we use expressions derived for time-dependent potentials, as expected, the deviations 
are systematically reduced when we increase level p. 

To show that the derived short-time effective actions can be used for calculating long-time 
transition amplitudes using the standard time discretization approach for path integrals, we 
display in Fig. |2] the convergence of the discretized long-time transition amplitude for the 



o 




A' 



FIG. 2. Convergence of numerical Monte Carlo results for the transition amplitude Ap(0,2;l,3) 
as a function of the number of time steps for the time-dependent harmonic oscillator ()6.2p with 
uj = 1, calculated with level p = 1,2,3 effective actions. The dashed lines give the fitted functions 
Ap + Bp/N'P + . . ., where the constant term Ap corresponds to the continuum-theory amplitude 
Ap{{), 2; 1, 3). The number of MC samples was N^c = 2 • 10^. MC error bars of around 2 x 10"^ 
for all values of p are shown in the graph. 
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Grosche-rescaled harmonic oscillator as a function of the number of time steps A^. As 
expected, we have obtained l/A^^ behavior, in accordance with our earlier conclusion on the 
A^-scaling of errors in this case. The numerical results presented in this graph are obtained 
using the Monte Carlo SPEEDUP [65] code, which was modified to include time-dependent 
effective actions implemented to high orders in the C programming language. 

The second exactly solvable model, which we have considered, is the forced harmonic 
oscillator |6lj . 

^FHo(2;,i) = 2'^^'^^ — xsini7t, (6.3) 

where Vt denotes the frequency of the external driving field. Fig. |3] presents numerical 
results for this model for = i7 = 1. The top plot gives deviations for the case of short-time 
transition amplitude, calculated analytically using level p effective action. We see again a 
perfect e-scaling of deviations. The bottom plot shows convergence of a long-time transition 



p = 2 — 



_ _ _ 

-m ■ l-— 



10 100 

A' 

FIG. 3. (left) Deviations |^p(2, 0; 2.001, e) — ^4(2, 0; 2.001, e)| as a function of propagation time 
e for the forced harmonic oscillator (j6.3p with a; = = 1, calculated analytically for p = 
1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20 from top to bottom. The dashed lines are proportional to e^'^ and 
e^^'^ , and demonstrate the perfect scaling of the corresponding level p = 1 and p = 20 results, 
(right) Convergence of numerical Monte Carlo results for the transition amplitude Ap(0, 0; 0.5, 1) 
as a function of the number of time steps N for p = 1,2,3. As before, dashed lines give the 
fitted functions Ap + Bp/N'^ + . . ., demonstrating the expected 1 /N'^ scaling of the deviations. The 
number of MC samples was Nmc = 2 • 10^. MC error bars of around 8 x 10^'' for all values of p 
are shown in the graph. 
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e 



amplitude as a function of the number of time steps A^, illustrating the expected 
behavior. 

We have numerically also considered the non-trivial case of a pure quartic oscillator 
V^q{x) = gx^/24: rescaled with the same Grosche factor = Vl + t^: 

^4 



Vc 



G,PQ 



[x,t) 



(6.4) 



24(l + t2)3 

This model is not exactly solvable, but the continuous transition amplitude can be de- 
termined numerically. To this end one could use the Grosche mapping, while relying on 
the previous numerical approach 
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providing exact transition amplitudes for the time- 



independent counterpart of the potential Vpq. Another possibility is based on numerically 
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FIG. 4. (left) Convergence of numerical Monte Carlo results for the transition amplitude 
Ap{0, 0; 0.5, 1) as a function of the number of time steps N for the pure quartic oscillator ()6.4p with 



the anharmonicity g = 0.1, rescaled by the Grosche factor ({t) = \/t^~+T, and calculated with 
level p = 1,2,3 effective actions. As before, dashed lines give the fitted functions Ap + Bp/N'^ + . . . . 
(right) In order to convincingly demonstrate the expected dominant 1 /N^ behavior of deviations, 
we plot I Ap(0, 0; 0.5, 1) — ^4(0, 0; 0.5, 1)| as a function of the number of time steps A^ for p = 1, 2, 3, 7. 
The dashed lines are fitted polynomials Ap + Bp/N^ + . . ., the same as on the top graph. The exact 
value of the amplitude is obtained as the constant term from fitting p = 7 results. The number of 
Monte Carlo samples was 

Nmc = 1.6- 10^3 for p = 7, in order to be able to resolve the exceedingly 
small deviations from the exact value of the amplitude. For p = 1, 2, 3 we used much smaller values 
of A^MC, typically 10*^ to 10^°. MC error bars of around 2 x 10"^° for p = 1, 2, 3 and 6 x 10"^^ for 
p = 7 are shown in both graphs. 
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obtained results from the modified SPEEDUP code, relying on the fitting of discretized am- 
plitudes to the e — )■ limit. As we can see from Fig. HI the numerical results exhibit a perfect 
scaling behavior for this non-trivial model as well. The power-law scaling of deviations in 
the calculated discretized transition amplitudes up to exceedingly small values fully verifies 
the presented analytic derivation of effective actions. 

As we have demonstrated during our analysis of several examples, the main advantage of 
the effective action approach is the calculation of transition amplitudes with high precision, 
which can be improved by using higher levels p. In practical applications, this increase 
in the precision is counterweighted by having to evaluate increasingly complex expression 
for the effective action. Fig. gives the comparison of the increase in the computational 
complexity of numerical codes for different levels p for the three models considered. As we 
can see, depending on the model, CPU time increases 5 to 10 times for p = 7, and clearly the 
benefit of 7 orders of magnitude increase in the precision of calculations far outweighs the 
CPU time increase. However, as p increases, the computational complexity will become too 
high. Since this strongly depends on the model, one has to perform a small-scale complexity 
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FIG. 5. Increase in the computational complexity of numerical simulations as a function of level 
p for the three models considered (Grosche-rescaled harmonic oscillator - G,HO; forced harmonic 
oscillator - FHO; Grosche-rescaled pure quartic oscillator - G,PQ). The complexity is measured by 
the relative increase in the CPU time Tcpu(p)/Tcpu(l) for calculation of transition amplitudes 
from Figs. [21 [3l and SI for each model correspondingly. The number of Monte Carlo samples was 
A^'mc = 2 • 10^ and number of time steps = 64. 
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study, similar to the one presented in Fig. [5l in order to find optimal value of level p to be 
used in numerical simulations. 

While the presented approach is excellently suited for applications where high-precision 
transition amplitudes are necessary, e.g. the calculation of partition functions, other higher- 



order schemes may be preferred 



split-operator method 
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'or solving different types of problems. For example, the 
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is ideally suited for studying both the real- and 
the imaginary-time evolution of various quantum systems. Furthermore, it has the advan- 
tage that it can be implemented in any chosen representation, which is appropriate for the 
quantum system, while the effective action approach relies on using the coordinate repre- 
sentation. 



VII. CONCLUSIONS 

We have presented an analytic procedure for deriving the short-time expansion of the 
propagator for a general one-particle non-relativistic quantum system in a time- dependent 
potential up to previously inaccessible high orders. The procedure is based on recursively 
solving both the forward and backward Schrodinger equation for the transition amplitude. 



Following an earlier approach for time-independent potentials j48| , we have derived recursion 
relations which allow an efficient analytic calculation of the effective potential to arbitrarily 
high orders in the propagation time e in the imaginary time. The analytically derived results 
are numerically verified by studying several simple models. 

The presented approach will be further expanded and generalized to many-body and 
multi-component systems in the forthcoming publication [l| . The next paper will also include 
generalization of the present imaginary-time approach to the real-time formalism, and its 
extensive numerical verification. 
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